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ABSTRACT 

Given the importance of shear flows for astrophysical gas dynamics, we study the 
evolution of the Kelvin-Hclmholtz instability (KHI) analytically and numerically. We 
derive the dispersion relation for the two-dimensional KHI including viscous dissi- 
pation. The resulting expression for the growth rate is then used to estimate the 
intrinsic viscosity of four numerical schemes depending on code-specific as well as on 
physical parameters. Our set of numerical schemes includes the Tree-SPH code VINE, 
an alternative SPH formulation developed by Price (2008) , and the finite- volume grid 
codes FLASH and PLUTO. In the first part, we explicitly demonstrate the effect of 
dissipation-inhibiting mechanisms such as the Balsara viscosity on the evolution of the 
KHI. With VINE, increasing density contrasts lead to a continuously increasing sup- 
pression of the KHI (with complete suppression from a contrast of 6:1 or higher). The 
alternative SPH formulation including an artificial thermal conductivity reproduces 
the analytically expected growth rates up to a density contrast of 10:1. The second 
part addresses the shear flow evolution with FLASH and PLUTO. Both codes result 
in a consistent non-viscous evolution (in the equal as well as in the different density 
case) in agreement with the analytical prediction. The viscous evolution studied with 
FLASH shows minor deviations from the analytical prediction. 

Key words: hydrodynamics - instabilities - methods:analytical - methods: numerical 
ISM:kinematics and dynamics 



1 INTRODUCTION 

Shear flows are an integral part of many astrophysical 
processes, from jets, the formation of cold streams, to 



outflows of protostars (Dekel et al. 2009 Agertz et al 



2009 Diemand et al. 2008 Walch et al. 20101, and cold 



gas clouds falling through the diffuse hot gas in dark 



matter halos (Bland-Hawthorn et al. 2007 Burkert et al 



2008). Jets and outflows of young stars can entrain ambient 
material, leading to mixing and possibly the generation of 
turbulence in e.g. molecular clouds (Burkert 2006; Banerjee 



et al. 2007 Gritschneder et al. 2009 Carroll et al. 20091, 



while the dynamical interaction of cold gas clouds with the 
background galactic halo medium can lead to gas stripping 
of e.g. dwarf spheroidals (e.g. Grcevich &i Putman 2009), 
and the disruption of high- velocity clouds ( Quilis & Moore 



vjunk@usm.uni-muenchen.de 



2001 Heitsch & Putman 20091. The KHI is believed to 



significantly influence the gas dynamics in all of these 
different scenarios. 

Moreover, viscous flows play a crucial role in e.g. gas 
accretion onto galactic disc s (|Das fc C hattopadhyay| |2008"l 
Park 2009 H einzeller et al.||2009[ ), as well as in dissipative 
processes like the turbulent cascade. Typically, the gas 
viscosity seems to be rather low in the interstellar medium, 
with typical flow Reynolds numbers of 10 5 . 
To describe these complex processes in detail, numerical 
schemes are applied to follow the hydrodynamical evolution. 
Numerous simulations use smoothed-particle hydrodynam- 



ics (SPH), (|Gingold fc Monaghan||1977| |Lucy||1977| |Benz 
|1990| |Monaghan| |1992| |2005| >, because its Lagrangian 
approach allows us to follow the evolution to high densities 
and small spatial scales. In combination with N-body 
codes, it is a perfect tool for cosmological simulations (e.g. 



Hernquist & Katz||1989| |Couchman et al.||1995; Springcl & 


Hernquist 2002 


Marri & White|2003 


Serna et al.|2003 


1 and 
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galaxy formation and evolution (Katz et al.| 1992 Evrard 



et al.|1994||Navarro et al.|1995| |Steinmetz fc Navarro|1999 



Thac ker fe Couchman| [20001 |Steinmetz fc Navarro| [2002 
Naab et al.||2006[ ). SPH describes the physical properties of 
a fluid by smoothing over a representative set of particles. 
However, this can lead to several problems. It can fail to 
correctly model sharp density gradients such as contact 
discontinuities, or velocity gradients occurring in e.g. shear 
flows (see Agertz et al. 20071, thus suppressing shear 



instabilities such as the KHI. 

An interesting problem to test the limitations of SPH as 
well as grid codes is the passage of a cold dense gas cloud 
moving through a hot and less dense ambient medium 



( Murray et al.||1993[ |Vietri et al.|p97l |Agertz et al.||2007 l 
Such a configuration would be typical for gas clouds raining 
onto galactic protodisks, for High- Velocity Clouds in the 
Milky Way and for cold HI clouds in the Galactic disk. 



Murray et al. ( 1993 1 demonstrated using a grid code that 



in the absence of thermal instabilities and/or gravity 
clouds moving through a diffuse gas should be disrupted by 
hydrodynamical shear flow instabilities within the time they 
need to travel through their own mass. Agertz et al. (2007) 
have shown that the KHI, and therefore the disintegration 
of such clouds is suppressed in SPH simulations. This 
problem, in particular the suppression of the KHI, has 
been subject to recent discussion in the literature. Several 
solutions have been proposed, e.g. Price (2008) discusses 



a mechanism, which involves a special diffusion term (see 
also |Wadsley et al. 12008]) 



Furthermore, Read et al. ( 2009 1 identify two effects occur 



ring in the SPH formalism, each one separately contributing 
to the instability suppression. The first problem is related 
to the leading order error in the momentum equation, which 
should decrease with increasing neighbor number. However, 
numerical instabilities prevent its decline. By introducing 
appropriate kernels, |Read et al.| ( |2009[ ) showed that this 
problem can be cured. The second problem arises due to the 
entropy conservation. Entropy conservation inhibits particle 
mixing and leads to a pressure discontinuity. This can be 
avoided by using a temperature weighted density following 



Ritchie & Thomas (2001). Recently, Abel (2010) has shown 



to solve this problem by evaluating the pressure force with 
respect to the local pressure. In contrast to standard SPH 
schemes this applies forces to particles only if there is a net 
force acting upon them. 

Another characteristic of SPH is the implementation of an 



artificial viscosity (AV) term (Monaghan & Gingold 1983), 
which is necessary in order to treat shock phenomena 
and to prevent particle interpenetration. AV can produce 
an artificial viscous dissipation in a flow corresponding 
to a decrease of the Reynolds-number and therefore a 



suppression of the KHI (Monaghan 20051. To confine this 
effect, a reduction of viscous dissipation was proposed 
by |Balsara| (|1995[) and improved by |Colagrossi| ( |2004[ ). 
Thacker et al. ( 2000 1 studied different AV-implementations 



in SPH and pointed out that the actual choice of the 
AV-implementation is the primary factor in determining 
code performance. An extension of SPH which includes 
physical fluid viscosities was d iscussed by e.g. |Takeda et al.| 
|T994|), |Flebbe et al.| (|l994| ) , |Espanol fc Revengaj (|2003[ ), 
Sijacki & Springel (20061 and Lanzafame et al. (20061. 



from a new class of hybrid schemes based on unstructured 
grids and combining the strengths of SPH and grid codes 



(Springel 2010). Some of the problems listed above might 
be solved with this type of implementation. 

In this paper we determine how accurate shear flows 
and the corresponding incompressible KHI are described in 
common numerical schemes. Therefore, in § [2j we analyti- 
cally derive the growth rates of the KHI including viscosity. 
In §[3] we briefly describe the numerical schemes and outline 
how the simulations have been analyzed. We then discuss 
our results. At first, we concentrate on the standard SPH 
implementation, which does not contain a physical viscosity 
but instead uses AV. However, as mentioned above, AV 
does influence the evolution of the flow. In § [4] we discuss 
the ability of two numerical SPH-schemes to model the 
incompressible KHI, namely the Tree-SPH method VINE 
(Wetzstein et al.||2009| |Nelson et al.||2009|), and the SPH 



code of Price (2008) 



By comparing to the derived analytical solution, we asses 
the effects of AV in VINE and estimate the intrinsic 



physical viscosity caused by AV (4.1l. We then study the 
development of the KHI for different density contrasts 



(4.2 I. We show that the instability is suppressed for density 
contrasts equal to or larger than 6:1. We also discuss the 
remedy suggested by Price (2008), hereafter P08. 
In Sj5]we then study the same problem with two grid codes, 
FLASH (|Fryxell et al.||2000|) and PLUTO JMignone et al. 



2007). As the intrinsic artificial viscosity is negligible in 



these schemes, we study the no n- viscous as well as the 



viscous evolution of the KHI for equal (5.1) as well as 



non-equal (5.2) density layers. We summarize our findings 
in m 



2 KHI - ANALYTICAL DESCRIPTION 

To derive the growth rate of the KHI in two dimensions 



including viscosity, we follow the analysis of Chandrasckhar 
( 1961 1 (for a related analysis see also Funada fc Joseph|2001 
and Kaiser et al. 20051. The fluid system is assumed to be 
viscous and incompressible. We use Cartesian coordinates in 
x and y, with two fluids at densities pi, p2, and velocities 
Ui, XJ% moving anti-parallel along the z-axis, separated by 
an interface layer at y — y s (see Fig. [I]). We neglect the 
effect of self-gravity. The hydrodynamical equations for such 
a system are then given by the continuity equation 

J^ + V-(pv) = 0, (1) 
and the momentum equation 



9v 

di 



+ (v ■ V) v 



= — Vp + pAv, 



(2) 



An alternative to conventional numerical schemes may arise 



with the flow density p, velocity v, the thermal pressure p 
and the kinematic viscosity v. 

2.1 Linear Perturbations 

We linearize equations [T] and [2] with the perturbations 

v -> v + (5v = (U(y) + u;w) (3) 
p -> p + Sp, (4) 
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Figure 1. Sketch of the initial conditions considered: Two fluid layers with constant densities pi and p2 flowing in opposite directions 
with uniform velocities Ui and U2- 



Sp. 



(5) 



u, w express the perturbation in the velocity, Sp and Sp in 
the density and pressure, respectively. This yields the system 
of linearized equations as 



pd t u + pUd x u + pwdyll = —d x Sp + u(p + 8p)d 2 U + 



pv(d x + dy)u, 



(6) 



pdtw + pUd x w 
d t 8p + Ud x Sp 
d t 5y B + U s d x Sy s 
5 x u + 5 y w 



-d y 8p + pv(d 2 x + dl)w, (7) 

-wd y p, (8) 

-vi(y.), (9) 

0. (10) 



Eqs. [6] and [7] represent the linearized Navier-Stokes equa- 
tions, where the density may change discontinuously at the 
interface positions denoted by y s . Eq.[8]is the linearized con- 
tinuity equation. In Eq. [9] the subscript s distinguishes the 
value of the quantity at y = y s (the interface layer). The 
last equation, Eq. [10] expresses the incompressibility of the 
fluid. With perturbations of the form 



u, w, Sp, Sp, 8y s ~ exp[i(k x x + fit)], 



(11) 



and assuming that the flow is aligned with the perturbation 
vector, i.e. k = k x , we arrive at 



D {p(n + kU)(Dw) - kp(DU)w} - pk 2 (n + kU)w 
iD{pvk 2 (Dw)} - iD {pv(D 3 w)} - 



D {kv(p + Sp)(D 2 U)} + ipuk 2 (D 2 w) - ipuk 4 

where D = d/dy. The term, ipvk 2 (D 2 w) in Eq. 
replaced with 

ipvk 2 (D 2 w) = ik 2 D(pv(Dw)) - ik 2 (Dw){D(pu)). 



12 



(12) 
can be 



(13) 



The boundary condition at y = y s is determined by an inte- 
gration over an infinitesimal element (y s — e to y s + e), for 
the limit e — > 0. Please note, that with Eq. [8] it follows for 
Sp, 



Sp- 



(14) 



After integration, the boundary condition becomes, 

A 3 {p(n + kU){Dw) - pk(DU)w} = 
ik 2 A s {vp(Dw)} - iA s {vp(D z w)} - 



kA a {vp{D 2 U)} -ikA a 



ik 2 A s {up(Dw)} — ik 2 lim 



(n + kU) 



(Dp)(D 2 U) + 



(Dw)D(up)dy (15) 



where A s is specifying the jump of any continuous quantity 
/ at y = y a , 



A s (/) = /(»=», +o) — /(»=», -o)- (16) 
For v = we retrieve the corresponding expression as given 



by Chandrasekhar (19611. 



2.2 Special case: constant velocities and densities 

To simplify the derivation of the growth rate n further, we 
consider the case of two fluid layers with constant densities 
pi and p2, and constant flow velocities U\ and U2 = —U\. 
In each region of constant pi, 2 and ?7i,2, Eq. |12| reduces to, 



[{n + kU ia )pi, 2 - 2ivk 2 \ (D 2 w) + iv(D 4 w) 
k 2 [(n + kU\;2) — ivk 2 ] w = 



(17) 



The layers are separated at y = y s = 0, and w/(n + kU) 
must be continuous at the interface. Also, w must be finite 
for y — ¥ 00, so that the solution of Eq. [T7]has the following 
form, 



w 
w 



A{n + kU l )e +ky (y < 0) 
A(n + kU 2 )e- hy (y>0). 



(18) 
(19) 



We assume that v\ = V2 = v (which is the case if we consider 
two media with the same viscous properties). Inserting this 
in Eq. |15| the characteristic equation yields, 



n + 2 



fc(Q2t/2 + OLlUl 



2 



n + 



(20) 
(21) 



k 2 (a 2 U% + a^l) - ik 3 u (a 2 U 2 + a%Ui) = 0. 

The parameters ot\, a.2 are defined by, 

P 1 P 2 
Qi = 1 , a 2 = 1 . 

pi + p2 pi + p2 

Solving for n, we get the expression for the mode of the 
linear KHI: 



(22) 



k(a2U2 + aiUi) 



ik v 



-k 2 aia 2 {Ui - U2) 2 



k*u 2 



(23) 



applying U2 = — Ui = U leads to 

ik 2 v 
IT 



k U (a.2 — ax) + 



-4fc 2 aia2(7 2 



4 



(24) 



The mode is exponentially growing/decaying with time, if 
the square root of n becomes imaginary, 

n = [k 2 U 2 (a.2 — evi)] + 
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physical parameters 


dimensionlcss 


in cgs units 


Box size 


2 


2 cm 


Mass 


4 


2780.81 g 


velocity 


0.387 


0.40 km/s 


time 


1 


9.8 -10- 6 s 



Table 1. Initial conditions in dimensionless units (first column) and in cgs units (second column). In the text we always refer to 
dimensionlcss units. 



vk 2 
2 



2 A: 4 



+ Ak 2 U 2 oaa2 



(25) 



The first term describes oscillations (which is not of interest 
for the growth), the second term the growth/decay, with a 
damping due to the viscosity. We use this formula for the 
comparison with our numerical studies for different density 
shearing layers. For equal density shearing layers pi = p2 — 
p, Eq. [25] leads to 



vk 2 
~2~ 



+ k 2 U 2 



(26) 



In §[4] and §[5] we use the velocity in direction of the pertur- 
bation, which in the above analysis refers to the y-direction 
and therefore, to the Uy-velocity component (w) when com- 
paring with simulations. The exponential term in Eq. [TT] 
(~ exp (i ■ n ■ t)) describes the time evolution of the KHI. In 
the following, we therefore compare In (v y ) with the analyt- 
ical expectation In (w) ~ i • n ■ t. 



3 KHI - NUMERICAL DESCRIPTION 

We use two independent numerical approaches - particle 
based and grid based - to follow the hydrodynamics of the 
system. In the following, all physical parameters are given 
in code units (see table 1 for conversion to physical units). 



3.1 SPH models - VINE & P08 



The parallel Tree -SPH code VINE ( |Wetzstein et al.||2009 
|Nelsonetal.|2009| has been successfully applied to a number 



of astrophysical problems on various scales (Naab et al 



2006 



Jesseit et al.||2007| |Gritschneder et al.||2009| |Walch et al.| 


2010 


Kotarba et al. 


2009 


1. In VINE the implementation 



of AV is based on the description by |Monaghan fe Gingold 



( 1983 1, and it includes the modifications by 



Lattanzio et al. 



( 19861. AV is not a real physical viscosity, but implemented 



to allow the treatment of shock phenomena. A viscous term, 

n 



n = -v 



eh 2 



(27) 



is added to the SPH momentum equations. The quantity 
e ~ 0.01 prevents a singularity if r — > 0, while h present the 
mean smoothing length between two particles. For v follows, 



■0- 



hv- 



eh? 



(28) 



p, and c are the mean density and the mean sound speed, 
respectively. The AV-parameter a controls the shear and 
the bulk viscosity, whereas the ft parameter regulates the 
shock-capturing mechanism. In the following we set a = 0.1, 
and j3 = 0.2 if not otherwise specified. AV reduces the 
Reynolds-number of the flow, resulting in the damping of the 
KHI p onaghan 2005). Balsara (1995) proposed a corrective 
term, improving the behavior of the AV in shear flows. Fur- 
ther improvements are discussed in Monaghan ( 2005[ ) and 
references therein. VINE can be run with and without the 
'Balsara- viscosity ' . 

To prevent the so-called 'artificial pairing' in SPH (e.g. 
Schuessler &i Schmitt||1981|), we implement a correction de- 



veloped by Thomas & Couchman (19921. Details can be 



found in Wetzstein et al. (2009) and 



Nelson et al. (20091. 



The SPH code presented in P08 uses a different implemen- 



tation of AV as explained in Morris ( 1997 1 to prevent the 



side effects of artificial dissipation. Additionally, a diffusion 
term called 'artificial thermal conductivity' is implemented 
(see § 4.2), which has been shown to prevent the KHI sup- 



pression in shear flows with large density contrasts (Price 



2008) 



3.2 Grid-based models - FLASH & PLUTO 

We choose the publicly available, MPI-parallel FLASH 
code version 2.5 ( |Fryxell et al.||2000[ ). FLASH is based on 
the block-structured AMR technique implemented in the 
PARAMESH library dMacNeice et al.||2000b. However, we 



do not make use of the AMR refinement technique, but use 
uniform grids throughout this paper. In FLASH'S hydrody- 
namic module the Navier-Stokes equations are solved using 
the piecewise parabolic method ( |Colella fc W oodward 1984] ), 
which incorporates a Riemann solver to compute fluxes be- 
tween individual cells. We use a Riemann tolerance value 
of 10~ 7 and a CFL of 0.5. Due to FLASH'S hydrodynamic 
scheme, the intrinsic numerical viscosity is reduced to a min- 
imum. This allows us to study the influence of a physical 
viscosity on the growth of the KHI. We therefore modify 
the hydrodynamical equations based on the FLASH module 
'diffuse' to explicitly include a viscous term, which scales 
with a given kinematic viscosity (see |5.1| and |5.2[ ). 
As an additional test, we apply the Godunov-type high res- 



olution shock capturing scheme PLUTO (Mignonc ct al. 
|2007[ ). It is a multiphysics, multialgorithm modular code, 
especially designed for the treatment of discontinuities. For 
the simulations described in this paper, we employ different 
Riemann-solvers and time-stepping methods on a uniform, 
static grid. 
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3.3 Initial conditions and analysis method 

Our numerical ICs are identical to the ones used for the 
derivation of the analytical growth rates (see Fig. [I] and 
table 1). To excite the instability, we apply a velocity per- 
turbation in y direction: 



v y — vo sin(fc • a;) ■ exp 



(29) 



where k is the wavenumber and vq is the perturbation 
amplitude of the y-velocity triggering the instability. 
The parameter controls how quickly the perturbation 
decreases with y (see discussion Appendix [A| . It is set 
to ctq = 0.1 if not otherwise specified. Initial pressure 
and density are set to po = 1 and po = 1, resulting in a 
sound speed of c 3i o = \/5/3 with an adiabatic exponent 
of 7 = 5/3. Since the analysis of s|2] is only valid for an 
incompressible fluid, the flow speed U must be subsonic. We 
chose U = 0.3 x c Sj o ~ 0.387, and the initial perturbation is 
Wo = 0.1 X U. We tested the assumption of incompressibility 
by calculating V • v, which vanishes for incompressible flows. 
This is satisfied in the linear regime, the primary focus of 
our work. The wavenumber k is equal to 4n/L, where L 
is the box length. The simulated box ranges from [—1, 1] 
in both directions. We use periodic boundary conditions. 
If not otherwise specified the AV parameters are set to 
a = 0.1 and /3 = 0.2. 

To analyze the SPH and grid simulations consistently, 
we bin the SPH particles on a 64 2 grid, using the cloud- 



in-cell method (Hockney & Eastwood 19881. For the grid 



codes, the same initial conditions are used. A resolution of 
512 2 is adopted during the calculation, but we rebin to a 64 2 
grid for the analysis. We measure the fastest-growing mode, 
which is the k — 4n/L mode of the velocity perturbation in 
y direction via a Fourier analysis. For more information see 
Appendix [B] 

We perform two sets of simulations with (i) equal density 
layers (see § 4.1 for SPH and § 5.1 for grid codes) and (ii) 



unequal density layers (see § 4.2 for SPH and § 5.2 for grid 
codes). In the latter case we assume pressure equilibrium. 
For SPH, we investigate the effects of equal mass and 
different mass particles (see 



4.2| 



4 SPH-SIMULATIONS OF THE KHI 

In the following, we model the evolution of the KHI in 
systems with pi = p2 ({ 4.1 1 and pi / p2 (j|4.2[). We ap- 



ply VINE, if not otherwise specified, and use the analytical 
growth rates (Eqs. 25 and 26 1 derived in § [5] to determine 
the effect of AV. 



4.1 Fluid layers with equal densities: 

In the case of pi = p2 we vary the following parameters: 
the resolution, which can be either enhanced by using more 
particles, or decreasing the smoothing length h, and the AV- 
parameters a and j3. We vary one parameter at a time, while 



context of AV we discuss the importance of the Balsara- 
viscosity. In Appendix [A] we also discuss the influence of 
different o"o, which determines the strength of the initial vy- 
perturbation (Eq. 291. 



• Dependence on resolution: 

According to the smoothing procedure in the SPH scheme, 
each particle requires a certain number of neighboring par- 
ticles for the calculation of its physical quantities. In VINE, 
these range from n neigh ,,nin to n neigh , m ax- The correspond- 
ing mean value of neighbors, n neig h, determines the smooth- 
ing length h. For a constant particle number, increasing 
n-neigh leads to a larger smoothing length, while at the same 
time the effective resolution is decreased. 
In Fig. [2] we show the time evolution of the u^-amplitude, 
which describes the growth of the KHI. For t ^ 0.2 the am- 
plitudes decrease since the SPH particles lose kinetic energy 
by moving along the y-direction into the area of the oppo- 
site stream (see Appendix [A| . Therefore we only consider 
t > 0.2 when fitting the growth rates of the KHI. The left 
panel of Fig.[2]shows the amplitude growth for n ne i g h = 20, 
30, and 40, respectively. (The commonly used value in two 
dimensions is n neig h = 30). All three cases appear to be sim- 
ilar. Thus, different n n ei g h do not have a substantial impact 
on the KHI-amplitude growth. 

The right panel of Fig. [2] shows the dependence on particle 
number, for the fiducial case of 512 2 (dotted line) and for 
an increased resolution of 1024 2 (solid line). The difference 
for the fitted viscosity is small (^ 1%). 

• Dependence of KHI on a, f3: 

In Fig. [3] Fig. [4] and Fig. [5] we show the KHI-evolution for 
different values of a and f3 without the Balsara-viscosity. 
Increasing the AV-parameter a or /3 results in a successive 
suppression of the KHI. Values of a > 2 and /3 > 1 lead to a 
decay of the initial perturbation. However, ft does not affect 
the growth as much as a. Therefore, we first concentrate on 
a as the operating term on the KHI. 

Can we assign an equivalent physical viscosity fsPH to the 
SPH scheme, i.e. can we determine how "viscous" the fluid 
described by SPH is intrinsically? To quantify its value, the 



the other ones are set to the fiducial values (see 3.1 1. In the 



analytical slope (Eq. 261, with the viscosity being the free 
parameter, is fitted to the simulated growing amplitudes. 
We show the best fits for a — 0.125 and a = 2 in the left 
panel of Fig. [5] for which we find the intrinsic viscosity of 
jy SPH = 0.07 and ^sph = 0.1. Here we assumed the time 
range of [0.2, 1], for which we determine the fits, to be well 
in the linear regime. 

In Fig.[6]we present the derived values of ^sph as a function 
of a. In summary, ^sph increases linearly with increasing 
plotted box size is from [—0.5,0.5] in both directions, the 
resolution is 512 2 . a, and the corresponding slope is 0.039. 
We also derive an offset of 0.065, which is the remaining 
intrinsic viscosity for a = 0. For each simulation, we also 
show the effective Re number of the flow (see Fig. |6j right 
y-axis), which was computed from Re = L ■ U/usph- The 
parameter L describes the characteristic scale of the pertur- 
bation, in our case the wavelength and U is the velocity of 
the flow. Clearly, the Reynolds-numbers we reach with our 
models are well below the commonly expected numbers for 
turbulent flows (Re > 10 5 ). 

The effective viscosity of the flow is also influenced by dif- 
ferent values of f3. Changing j3 by a factor of two (e.g. from 
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Figure 3. Time evolution of the KHI using VINE for increasing AV parameter a (top to bottom) and constant = 2 The panels show 
the central region of each simulation box, ranging from [—0.5,0.5]. The upper layer (grey area) is moving to the left, the lower layer 
(black area) to the right. Noticeable damping occurs for a > 0.125 (see left panel of Fig. pi. 



P = 0.5 to /3 = 1) results in an increase in effective viscosity 
by a factor of 0.01 (see right panel of Fig. [5|. 

• Dependence on the Balsara-viscosity: 
We showed that AV leads to artificial viscous dissipation, 
resulting in the damping of the KHI. To prevent this, we 
use the Balsara-viscosity, see also section |3.1| In Fig. [7] we 
show the corresponding amplitudes for three examples of 
AVs: (a = 0.1, /3 = 0.2), (a = 1, P = 2) and (a = 2, P = 2). 
Clearly, the Balsara viscosity reduces the damping of the 
KHI, rendering ^sph almost independently of a and /3 (see 
also Fig. [6}. 



4.2 Fluid layers with variable densities: 

While the previously addressed case of equal densities 
helped us to understand the detailed evolution of the KHI 
as modeled with SPH, the astrophysically more interesting 
case are shear flows with different densities. The resolution 
of the diffuse region is lower by a factor of \/ DC, where DC 
is the ratio of the densities in dense and diffuse medium 
(e.g. DC = 10 corresponds to a density contrast of 10 : 1). 
We return to our standard set of parameters, in which case 
a = 0.1 and /3 = 0.2. For these low AV parameters we do 
not need the Balsara-viscosity (see 4.1 1. (Nonetheless, we 
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Figure 4. Like Fig. [3] but for increasing values of the AV-parameter j3 (a 
panel of Fig. [si. 



0.1). A noticeable damping occurs for of ft > 1 (see right 





Figure 5. Left panel: Time evolution of the VINE Wy-amplitude for different values of the AV-parameter a, where f3 has been fixed to 
/3 = 2. The thick dashed-dotted lines correspond to the analytical fit, shown for a = 0.125 and a = 2 (which corresponds to iajph = 0.07 
and i^sph = 0.1). Right panel: Like before, but for different values of the AV-parameter /3, where a has been fixed to a = 0.1. 



did run test simulations with the Balsara switch, which we 
found to confirm our former finding, since the growth of 
the KHI was not affected). In the following, we (i) analyze 
the growth of the KHI for different values of DC (with 
equal mass particles) and address the problem of KHI 
suppression, while in (ii) we test the influence of equal mass 
or spatial resolution. 

(i) KHI GROWTH AS A FUNCTION OF DC: 
We show the KHI evolution for increasing DC in Fig. [8] 
For DC ^ 6 the KHI does not develop anymore. This 
SPH problem of KHI suppression has been studied in great 



detai l (e.g. |Agertz et al.||2007| |Price||2008| |Wadsley et~ 
2008}|Read et al.|2009||Abel|2010| ). SPH particles located at 



the interface have neighbors at both sides of the boundary 
(i.e. from the dense- and less dense region). Therefore, the 
density at the boundary is smoothed during the evolution. 
However, the corresponding entropy (or, depending on 
the specific code, the thermal energy) is artificially fixed 
in these (isothermal) setups which results in an artificial 
contribution to the SPH pressure force term, due to which 
the two layers are driven apart. One possible solution is to 



either adjust the density (Ritchie & Thomas 2001 Read 



et al. 20091, or to smooth the entropy (thermal energy) 
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Figure 6. Derived physical viscosities (jajph) corresponding to different AV parameters a with (open red points) and without (filled 
black points) Balsara-viscosity. We also show the corresponding effective Re-numbers. 



a = o.i /S= 0.2 

a = 1.0 0= 2.0 




-2 h d 
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time 

Figure 7. Time evolution of the VINE Dy-amplitude for different values of the AV-parameters a and fi, where the Balsara-viscosity has 
been used. The damping of the amplitudes is completely prohibited by the Balsara switch. 



( Price|2008||Wadsley et al.|2008||Abel|2010| . 

A remedy has been discussed by |Price| ( |2008| , who proposed 
to add a diffusion term, which is called artificial thermal 
conductivity (ATC), to adjust the thermal energy. (For a 



detailed study of ATC see Price 2008). With this method, 



the KHI should develop according to the test cases of P08. 
In Fig. [9] we test whether the P08 approach is indeed in 
agreement with our analytical prediction. Note that P08 has 
a method implemented to account for the artificial viscous 
dissipation caused by AV (similar to the Balsara-viscosity). 
Thus, the viscous effects of AV are strongly reduced. For 
DC = 10 and using 512 2 particles in the dense layer 
we indeed find good agreement between measured and 
analytical growth rates. If the standard SPH scheme is 
used, a correction term like ATC has to be included to 



obtain a KHI in shear flows with different densities, which 
is consistent with the analytical prediction. 

(ii) KHI GROWTH USING EQUAL AND DIFFERENT PAR- 
TICLE MASSES: 

First, we investigate the development of the KHI for the 
standard SPH case of equal mass resolution throughout 
the computational domain, and therefore fewer particles 
in the low density fluid layer (see top panel of Fig. 



10 for 



DC = 10, where the dense medium is resolved with 512 2 
particles). This results in a varying spatial resolution, due 
to the fact that SPH derives the hydrodynamic quantities 
within a smoothing length h set by a fixed number of 
nearest neighbors. This construct - as has been discussed 
in detail earlier in e.g. Agertz et al. 2007 - specifically 
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Figure 8. Like Fig. |10| top panel, but for different density contrasts. From top to bottom we show DC=2, 3, 6. For DC ^ 6 the KHI 
does not develop anymore. 



PRICE for 10:1 
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Figure 9. Time evolution of the KHI modeled with P08 for the DC = 10. The dashed-dotted line corresponds to the analytical 
prediction, Eq. |30| which is in good agreement with the simulation. 



lowers the Reynolds-number of the shear flow across density 
discontinuities, thus affecting the evolution of the KHI. 
As can be seen in the top panel of Fig. |10| the KHI is 
completely suppressed. 

Second, we test the case of equal spatial resolution in 
both fluid layers, and therefore unequal particle masses 
within the computational domain (Fig. 10 lower panel). 



Again, we find the KHI to grow too slowly with respect to 
the analytical estimate. However, the suppression is less 
effectively in the latter case. 



5 GRID-SIMULATIONS OF THE KHI 

For comparison to the SPH treatment of Kelvin-Helmholtz 
instabilities, we study an identical setup of fluid layers with 
the grid-based codes FLASH and PLUTO (see jpL2}. We 



reuse the previously specified initial conditions with a grid 
resolution of 512 2 cells in the standard case. For FLASH, we 
additionally include physical viscosity of various strenghth 
in some of the simulations (see j ]3.2[ ). Note, that for the 
following examples we use ao — 1 if not otherwise specified, 
which does not affect the growth of the amplitudes in the 
linear regime (for further information see discussion in the 
Appendix [Ah. 
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Figure 10. KHI with VINE for DC = 10. Top: Case of equal particle masses. Bottom: Case of unequal particle masses and therefore 
equal particle numbers in both layers. The KHI is suppressed in all cases. 
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Figure 11. Evolution of KHI amplitudes for equal density layers. Left panel: Non-viscous evolution for FLASH (solid line), and PLUTO 
(dotted line). Additionally, we show the example with VINE (dashed line), where the AV has been set to zero (a = = 0). Right panel: 
Viscous KHI evolution using FLASH. The thick dashed-dotted lines correspond to the analytical prediction, Eq. |26| 



5.1 Fluid layers with equal densities 



5.1.2 Viscous evolution 



5.1.1 Non-viscous evolution 



The left panel of Fig. [TT] shows the non-viscous KHI- 
evolution, using FLASH (solid line), PLUTO (dotted line), 
and for comparison VINE (dashed line). In the VINE ex- 
ample, the AV has been set to zero (a = ft = 0). The 
expected analytical growth (Eq. 26 1 reduces with v = 
to n ~ k ■ U = 2.43 (indicated by the thick dashed dotted 
line). The FLASH and PLUTO amplitudes develop in a sim- 
ilar pattern and are almost undistinguishable. Their fitted 
slopes within the linear regime (which lies roughly between 
t = 0.3 - 0.6) results in n flt = 2.49. FLASH and PLUTO 
show a consistent growth in agreement with the analyti- 
cal prediction. VINE on the other hand exhibits a slightly 
slower growth. This deviation is due to the intrinsic viscosity 
(Vint = 0.065) that was estimated in 4.1 



The right panel of Fig. |ll| shows the viscous KHI-amplitudes 
using FLASH. The corresponding analytical predictions 
(Eq. 26 1 are shown by the thick dashed-dotted lines for 
the examples with v = 0.00003 and v — 0.03. To quan- 
tify the growth of the KHI in the FLASH simulations, we 
again fit the slopes of the KHI-amplitude in the linear regime 
(between t — 0.3 — 0.6). The result along with the corre- 
sponding error is plotted in Fig. |12| For small viscosities 
(y < 0.003), we find the growth rates of the KHI in FLASH 
to be in good agreement with the analytical prediction. In 
this viscosity range, the dominant term in the analytical 
prediction (Eq. 26 1 is ~ kU. Therefore, any influence of v 
is marginal, and the amplitudes do not change considerably. 
FLASH treats the fluid as if v w 0. 

However, with increasing viscosity, the amplitudes should be 
damped. This behavior is in fact visible in the right panel of 
Fig. [TT] (as well as in Fig. [ljf . The growth rates of the KHI 
agree very well with the analytical prediction. 
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Figure 12. Comparison of the analytical expectation and the models for DC=1 (diamond shaped symbols) and DC = 10 (square 
symbols). The slopes derived for FLASH correspond to the analytical fits. The lines represent the analytic prediction, for DC=1 (solid 
line, see Eq. |26|l and DC = 10 (dashed line, see Eq. |25|l. 




Figure 13. Time evolution of the KHI density in a simulation with nu=0 and DC = 10 for FLASH (top row) and PLUTO (bottom 
row). The plotted box size is from [—1, 1] in both directions, the resolution is 512 2 . The KHI develops, which is in contrast to the example 
simulated with VINE. 



5.2 Fluid layers with different densities 

5.2.1 Non-viscous evolution 

Finally, we investigate a density contrast of 10 : 1, similar 



to the example studied with VINE (see § 4.2 1. Fig. 13 shows 



the non-viscous evolution of the KHI for the DC = 10 case 
(upper line for FLASH, bottom line for PLUTO). It can be 
seen that for both codes the interface layer starts to roll-up 
and the instability is developed. This is in disagreement with 
the previously discussed case using SPH, where the KHI is 
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Figure 14. The same as in Fig.[TT]but for a DC = 10. Left panel: Non-viscous evolution for FLASH (solid line), PLUTO (dotted line) 
and the high-resolution (1024 2 ) amplitude for FLASH (dashed line). Right panel: Non-viscous evolution using PLUTO, with different 
solvers, see text for more details. 
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Figure 15. Viscous evolution using FLASH. The thick dashed-dotted lines correspond to the analytical prediction, Eq. |25| 



completely suppressed for DC> 6 (see |4.2[ ). 
The left panel of Fig. [14] presents the corresponding am- 
plitudes for FLASH (solid line) and PLUTO (dotted line) 
compared to the analytical prediction (thick dashed-dotted 
line), which in this case reduces to 



= ±i-s/4k 2 U 2 aiot2- 



(30) 



For FLASH we show two different resolutions (512 2 and 
1024 2 ). The amplitudes resulting in the case of low and 
high resolution are effectively indistinguishable. This is an 
important result, as it demonstrates that small scale per- 
turbations, which arise due to numerical noise and which 
could violate the linear analysis (as we then might follow 
the growth of higher order modes rather than the initial 
perturbation) are not important. Therefore, we have shown 
that our simulations are converged as we would otherwise 
expect the growth of the KHI to be slightly dependent on 
the grid resolution (see e.g. the recent findings of |Robert^| 
son et al. (20101, who had to smooth the density gradient 



between the two fluid layers in order to achieve convergence 
in terms of grid resolution). Moreover, both FLASH and 
PLUTO evolve similarly. For all three examples the slope of 



the amplitude evolution can be approximated to 1.4, which 
is in good agreement with the analytical expectation. Note 
that we do not show the comparison with the VINE ampli- 
tude since the KHI does not evolve for DC — 10 (see 4.2 1. 



Many grid codes offer a variety of hydrodynamical solvers. 
We therefore tested the influence of different numerical 
schemes on the growth of the KHI using PLUTO (see 
right panel of Fig. 



14 



using 

We show three different examples; 



'simOOO' is a Lax-Friedrichs scheme together with a sec- 
ond order Runge-Kutta solver (tvdlf); 'simOOl' implements 
a two-shock Riemann solver with linear reconstruction em- 
bedded in a second order Runge-Kutta scheme; 'sim002' also 
implements a two-shock Riemann solver, but with parabolic 
reconstruction, and embedded in a third order Runge-Kutta 
scheme. Both, 'simOOl' and 'sim002' show a similar growth 
of the KHI in agreement with the analytical prediction (see 
Fig. 14 top right panel). The more diffusive scheme used 
in 'simOOO' causes a small delay in the growth of the KHI, 
but results in a similar slope within the linear regime (up to 
t = 0.6). 
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5.2.2 Viscous evolution 

Fig. [15] shows the viscous KHI-amplitudes using FLASH, 
which are increasingly suppressed with v. The correspond- 
ing analytical prediction (Eq. 25 1 is shown for v — 0.0003, 
and v — 0.03 (thick dashed-dotted lines). For v < 0.03 the 
simulated growth rate is slightly enhanced by a factor of 
~ 0.12 as compared to the analytical prediction (see also 
Fig. 121. However, for higher viscosities (v ^ 0.03) we find 
good agreement between simulation and analytical predic- 
tion. 



6 CONCLUSIONS 

We have studied the Kelvin-Helmholtz instability applying 
different numerical schemes. We use two methods for our 



between the analytically predicted amplitude evolution and 
the simulation oflPricel (120081) for DC = 10. 



et al. 



Price 



SPH models, namely the Tree -SPH code VINE fWetzstein 
2009||Nelson et al.||2009[ ), and the code developed by 



(20081). The grid based simulations of the KHI rely on 
FLASH ( |Fryxell et al.||2000[ ), while as a t est for the non- 
viscous evolution we also apply PLUTO (Mignone et al 



2007) 



We first extended the analytical prescription of the KHI by 
|Chand rasckhar (1961 1 to include a constant viscosity. With 
this improvement we were able to measure the intrinsic vis- 
cosity of our subsequently performed numerical simulations. 
We test both SPH as well as grid codes with this method. 
We then concentrated on the KHI-evolution with SPH. We 
performed a resolution study to measure the dependence 
of the KHI growth on the mean number of SPH neighbors 
(usph) and the total number of particles, respectively. We 
found that our simulations were well resolved and that a 
different number of usph did not significantly influence the 
KHI growth rate. 

In case of equal density shearing layers we then measured 
the intrinsic viscosity in VINE by evaluating our simulations 
against the analytical prediction in the linear regime. With- 
out using the Balsara viscosity the AV parameters a and 
P effectively lead to a damping of the KHI. The commonly 
suggested and used settings of a = 1, and /3 — 2 result in a 
strong suppression of the KHI. More quantitatively, we de- 
rive values of 0.065 < ^sph < 0.1 for < a < 1. Different 
values of j3 do not have a strong impact on ^sph- By in- 
troducing the Balsara- viscosity the dissipative effects of the 
AV can be reduced significantly, effectively rendering the re- 
sults to be independent of a and /3. However, the constant 
floor viscosity of ^sph = 0.065 prevails. Furthermore for a 
given a, we estimated the effective Reynolds-number (Re) of 
the flow. For the minimum SPH viscosity of mtspH = 0.065 
we derive a maximum Reynolds number of 12. This is very 
small compared to typical Reynolds numbers of real turbu- 
lent flows (Re > 10 ). For different density shearing layers 
we confirmed the results discussed inlAgertz et al. (20071, i.e. 



the KHI is completely suppressed for shear flows with differ- 
ent densities (in the case of VINE for DC ^ 6). Here, using 
the Balsara switch does not solve the problem. This indicates 
that other changes to the SPH formalism are required in or- 
der to correctly model shearing layers of different densities. 



To demonstrate this we applied the solution of Price ( 2008 1 
to our initial conditions for DC=10. In this case the KHI 
was suppressed in VINE. However, we found good agreement 



The second part of this paper addresses the non-viscous- 
and viscous KHI evolution using grid codes. In the case 
of equal density shearing layers, we found the non-viscous 
growth rates for shear flows with FLASH and PLUTO to 
be in good agreement with the analytical prediction. In the 
viscous case, the FLASH-amplitudes show only a minor de- 
pendency on the viscosity if v < 0.03. Increasing the viscos- 
ity leads to a damped evolution, with the simulated growth 
coinciding with the analytical prediction. 
For non-viscous shear flows (with a density contrast of 
DC = 10) the KHI does develop for FLASH and PLUTO in 
agreement with the analytical prediction. In the viscous case 
FLASH (also analyzed with DC = 10) slightly overpredicts 
the corresponding growth rates for v < 0.03 by a constant 
factor of ~ 0.12. 

The comparison between VINE, FLASH and PLUTO in the 
equal density case, where AV = and v = 0, demonstrated 
that VINE does have an intrinsic viscosity (which we esti- 
mated to ~ 0.065). 



ACKNOWLEDGMENTS 

I would like to thank Volker Springel for his useful sugges- 
tions, and Thorsten Naab for his support. Many thanks also 
to Oscar Agertz for helpful discussions, as well as Eva Ntor- 
mousi. This research was supported by the DFG priority 
program SPP 1177 and by the DFG cluster of excellence 
'Origin and Structure of the Universe'. Part of the simula- 
tions were run on the local SGI ALTIX 3700 Bx2 which was 
also partly funded by this cluster of excellence. FLASH was 
developed by the DOE-supported ASC/ Alliance Center for 
Astrophysical Thermonuclear Flashes at the University of 
Chicago. S.Walch gratefully acknowledges the support of the 
EC-funded Marie Curie Research Training Network Con- 
stellation (MRTN-CT-2006-035890). M. Wetzstein grate- 
fully acknowledges support from NSF grant 0707731. 



REFERENCES 

Abel T., 2010, ArXiv e-prints 

Agertz O., Moore B., Stadel J., Potter D., Miniati F., Read 
J., Mayer L., Gawryszczak A., Kravtsov A., Nordlund 
A., Pearce F., Quilis V., Rudd D., Springel V., Stone J., 
Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MN- 
RAS, 380, 963 

Agertz O., Teyssier R., Moore B., 2009, MNRAS, 397, L64 
Balsara D. S., 1995, Journal of Computational Physics, 121, 
357 

Banerjee R., Klessen R. S., Fendt C, 2007, ApJ, 668, 1028 
Benz W., 1990, in Buchler J. R., ed., Numerical Modelling 

of Nonlinear Stellar Pulsations Problems and Prospects 

Smooth Particle Hydrodynamics - a Review, pp 269 — V 
Bland-Hawthorn J., Sutherland R., Agertz O., Moore B., 

2007, ApJL, 670, L109 
Burkert A., 2006, Comptes Rendus Physique, 7, 433 
Burkert A., Naab T., Johansson P. H., Jesseit R., 2008, 

ApJ, 685, 897 



14 Junk et al. 



Carroll J. J., Frank A., Blackman E. G., Cunningham A. J., 
Quillen A. O, 2009, ApJ, 695, 1376 

Chandrasekhar S., 1961, Hydro dynamic and Hydromag- 
netic Stability. Dover Publications 

Colagrossi A., 2004, Dottorato di Ricerca in Meccanica 
Teorica ed Applicata XVI CICLO, A meshless Lagrangian 
method for free-surface and interface flows with fragmen- 
tation, Phd Thesis. Universita di Roma, La Sapienza 

Colella P., Woodward P. R., 1984, Journal of Computa- 
tional Physics, 54, 174 

Couchman H. M. P., Thomas P. A., Pearce F. R., 1995, 
ApJ, 452, 797 

Das S., Chattopadhyay I., 2008, New Astronomy, 13, 549 
Dekel A., Birnboim Y., Engel C, Freundlich J., Goerdt 

T., Mumcuoglu M., Neistein E., Pichon C, Teyssier R., 

Zinger E., 2009, Nature, 457, 451 
Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., 

Potter D., Stadel J., 2008, Nature, 454, 735 
Espanol P., Revenga M., 2003, Phys. Rev. E, 67, 026705 
Evrard A. E., Summers F. J., Davis M., 1994, ApJ, 422, 

11 

Flebbe O., Muenzel S., Herold H., Riffert H., Ruder H., 

1994, ApJ, 431, 754 
Fryxell B., Olson K., Ricker P., Timmes F. X., Zingale M., 

Lamb D. Q., MacNeice P., Rosner R., Truran J. W., Tufo 

H., 2000, ApJS, 131, 273 
Funada T., Joseph D. D., 2001, Journal of Fluid Mechanics, 

445, 263 

Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375 
Gritschneder M., Naab T., Burkert A., Walch S., Heitsch 

F., Wetzstein M., 2009, MNRAS, 393, 21 
Gritschneder M., Naab T., Walch S., Burkert A., Heitsch 

F., 2009, ApJL, 694, L26 
Heinzeller D., Duschl W. J., Mineshige S., 2009, MNRAS, 

397, 890 

Heitsch F., Putman M. E., 2009, ApJ, 698, 1485 

Hernquist L., Katz N., 1989, ApJS, 70, 419 

Hockney R. W., Eastwood J. W., 1988, Computer simula- 
tion using particles. Institute of Physics Publishing 

Jesseit R., Naab T., Peletier R. F., Burkert A., 2007, MN- 
RAS, 376, 997 

Kaiser C. R., Pavlovski G., Pope E. C. D., Fangohr H., 

2005, MNRAS, 359, 493 
Katz N., Hernquist L., Weinberg D. H., 1992, ApJL, 399, 

L109 

Kotarba H., Lesch H., Dolag K., Naab T., Johansson P. H., 

Stasyszyn F. A., 2009, MNRAS, 397, 733 
Lanzafame C, Belvedere G., Molteni D., 2006, AAP, 453, 

1027 

Lattanzio J., Monaghan J., Pongracic H., Schwartz M., 
1986, SIAM Journal on Scientific and Statistical Com- 
puting, 7, 591 

Lucy L. B., 1977, AJ, 82, 1013 

MacNeice P., Olson K. M., Mobarry C, de Fainchtein R., 
Packer C, 2000, Computer Physics Communications, 126, 
330 

Marri S., White S. D. M., 2003, MNRAS, 345, 561 
Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu 

O., Zanni C, Ferrari A., 2007, ApJS, 170, 228 
Monaghan J. J., 1992, ARA& A, 30, 543 
Monaghan J. J., 2005, Reports of Progress in Physics, 68, 

1703 



Monaghan J. J., Gingold R. A., 1983, Journal of Compu- 
tational Physics, 52, 374 
Morris J., 1997, Journal of Computational Physics, 136, 41 
Murray S. D., White S. D. M., Blondin J. M., Lin D. N. C, 

1993, ApJ, 407, 588 
Naab T., Jesseit R., Burkert A., 2006, MNRAS, 372, 839 
Navarro J. F., Frenk C. S., White S. D. M., 1995, MNRAS, 
275, 56 

Nelson A. F., Wetzstein M., Naab T., 2009, ApJS, 184, 326 
Park M., 2009, ApJ, 706, 637 

Price D. J., 2008, Journal of Computational Physics, 227, 
10040 

Quilis V., Moore B., 2001, ApJL, 555, L95 
Read J. I., Hayfield T., Agertz O., 2009, ArXiv 0906.0774 
Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743 
Robertson B. E., Kravtsov A. V., Gnedin N. Y., Abel T., 

Rudd D. H., 2010, MNRAS, 401, 2463 
Schuessler I., Schmitt D., 1981, AAP, 97, 373 
Serna A., Dommguez-Tenreiro R., Saiz A., 2003, ApJ, 597, 

878 

Sijacki D., Springel V., 2006, MNRAS, 371, 1025 
Springel V., 2010, MNRAS, 401, 791 
Springel V., Hernquist L., 2002, MNRAS, 333, 649 
Steinmetz M., Navarro J. F., 1999, ApJ, 513, 555 
Steinmetz M., Navarro J. F., 2002, New Astronomy, 7, 155 
Takeda H., Miyama S. M., Sekiya M., 1994, Progress of 

Theoretical Physics, 92, 939 
Thacker R. J., Couchman H. M. P., 2000, ApJ, 545, 728 
Thacker R. J., Tittley E. R., Pearce F. R., Couchman 

H. M. P., Thomas P. A., 2000, MNRAS, 319, 619 
Thomas P. A., Couchman H. M. P., 1992, MNRAS, 257, 

11 

Vietri M., Ferrara A., Miniati F., 1997, ApJ, 483, 262 
Wadsley J. W., Veeravalli C, Couchman H. M. P., 2008, 

MNRAS, 387, 427 
Walch S., Naab T., Whitworth A., Burkert A., 

Gritschneder M., 2010, MNRAS, 402, 2253 
Wetzstein M., Nelson A. F., Naab T., Burkert A., 2009, 

ApJS, 184, 298 



APPENDIX A: DEPENDENCE OF 
KHI- AMPLITUDES ON a 

Dependence of KHI on ctq: 

This parameter determines the strength of the initial v y - 
perturbation (Eq. 291. In Fig. [2] we show the time evolution 
of the vy-amplitude, which describes the growth of the KHI. 
For t ^ 0.2 the amplitudes decrease since the SPH particles 
lose kinetic energy by moving along the y-direction into the 
area of the opposite stream. If the magnitude of the initial 
perturbation is low (i.e. small <to), then the decrease in the 
amplitude is stronger than for e.g. (To = 1, where the initial 
perturbation is large and the decrease less prominent. But 
independently of the value of oo the subsequent growth of 
the instability is similar, and we obtain comparable results 
neglecting the decreasing initial part. Fig. |A1| shows the de- 
pendency of the KHI-amplitudes using different values of 
cro, for VINE (left side) and FLASH (right side). For this 
example we use equal density layers, where for FLASH a 
viscosity of v = 0.3 has been taken. Clearly visible is the 
initial drop caused by a low value of oq. This is the case for 
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Figure Al. Variation of KHI-amplitudc in the case of equal density layers using VINE (left side) and FLASH (right side) for different 
values of erg. For Flash the viscosity has been set to v = 0.3 



both codes, and arises due to the transformation of energy 
to build up the KHI. The fitted slopes do not vary much 
with (To. To extract the slopes, we concentrate on the time 
evolution after this initial drop. 



APPENDIX B: MEASURING THE 
KHI-AMPLITUDES 

To measure the amplitude growth of the KHI, we apply a 
Fourier- Transformation (FT) to the Vy-velocity component 
of the grid points. The FT allows to select the desired modes 
reducing the numerical noise. 

The region of our focus, x = [—0.5,0.5] and y = [—0.5,0.5] 



contains one mode of the u^-perturbation (Eq. 29 I triggering 
the instability, see Fig. |B1[ The shaded regions comprise the 
particles subject to the FT (at x = —0.5 and x = 0.5). The 
maximum of the FT gives the dominant mode k and its 
corresponding velocity amplitude, which we compare with 
the analytical model. 
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Figure Bl. Method to measure the KHI-amplitudes: The Dy-velocity of the particles within the shaded region are subject to the 
Fourier- Transformation. The maximum of the Back- Transformation gives the maximal amplitude. 



